function [betas_tightness, betas_f] = reg_hpxr_real_nom(gesol,Phi_NxYz,prob_NxYz,Y_grid,FLAG_RUN_PARALLEL)
% slope coefficients for the following regression:
% hpxr(n,t+12) = a + b*X(t) + e(t)
% hpxr(n,t+12) = p(n-1,t+12) - p(n,t) + p(1,t)
% X(t) is either tightness(t) or f(t)
% maturity is n=2,3,4,5 years

    if nargin<5; FLAG_RUN_PARALLEL = false; end

    Phi_NxYz_12mo = mpower2(Phi_NxYz,12);
    
    betas_tightness = zeros(1,4);
    betas_f = zeros(1,4);

    if FLAG_RUN_PARALLEL
        
        % prepare inputs
        lhs0_list = cell(4,1);
        lhst_list = cell(4,1);
        rhs0a_list = cell(4,1);
        rhs0b_list = cell(4,1);
        
        inflation_process = gesol.inflation_process;
        
        for n = 2:5
            
            idx_buy = find(gesol.nom_bonds.maturities==12*n);
            idx_sell = find(gesol.nom_bonds.maturities==12*(n-1));
            idx_1y = find(gesol.nom_bonds.maturities==12);

            lhs0.fun_Nxz = -log(gesol.nom_bonds.P_A(:,:,:,idx_buy)) + log(gesol.nom_bonds.P_A(:,:,:,idx_1y));
            lhs0.slope_Y = gesol.nom_bonds.P_C_a(idx_buy) - gesol.nom_bonds.P_C_a(idx_1y);
            lhs0.const = gesol.nom_bonds.P_B_a(idx_buy) - gesol.nom_bonds.P_B_a(idx_1y);
            lhs0.slope_X = gesol.nom_bonds.P_B_b(idx_buy) - gesol.nom_bonds.P_B_c(idx_1y);
            lhs0.slope_eps = gesol.nom_bonds.P_B_c(idx_buy) - gesol.nom_bonds.P_B_c(idx_1y);
            
            lhs0_list{n-1} = lhs0;

            lhst.fun_Nxz = log(gesol.nom_bonds.P_A(:,:,:,idx_sell));
            lhst.slope_Y = -gesol.nom_bonds.P_C_a(idx_sell);
            lhst.const = -gesol.nom_bonds.P_C_a(idx_sell);
            lhst.slope_X = -gesol.nom_bonds.P_B_b(idx_sell);
            lhst.slope_eps = -gesol.nom_bonds.P_B_c(idx_sell);
            
            lhst_list{n-1} = lhst;

            % rhs = tightness
            rhs0.fun_Nxz = gesol.theta;
            rhs0.const = 0;
            rhs0.slope_Y = 0;
            rhs0.slope_X = 0;
            rhs0.slope_eps = 0;
            
            rhs0a_list{n-1} = rhs0;
            
            % rhs = f
            rhs0.fun_Nxz = gesol.f;
            
            rhs0b_list{n-1} = rhs0;
            
        end
        
        parfor n = 2:5

            [a,betas_tightness(n-1),R2] = generic_nominal_reg(inflation_process, lhs0_list{n-1}, 12, lhst_list{n-1}, rhs0a_list{n-1}, prob_NxYz, Phi_NxYz_12mo, Y_grid);

            [a,betas_f(n-1),R2] = generic_nominal_reg(inflation_process, lhs0_list{n-1}, 12, lhst_list{n-1}, rhs0b_list{n-1}, prob_NxYz, Phi_NxYz_12mo, Y_grid);

        end
        
    else
    
        for n = 2:5

            idx_buy = find(gesol.nom_bonds.maturities==12*n);
            idx_sell = find(gesol.nom_bonds.maturities==12*(n-1));
            idx_1y = find(gesol.nom_bonds.maturities==12);

            lhs0.fun_Nxz = -log(gesol.nom_bonds.P_A(:,:,:,idx_buy)) + log(gesol.nom_bonds.P_A(:,:,:,idx_1y));
            lhs0.slope_Y = gesol.nom_bonds.P_C_a(idx_buy) - gesol.nom_bonds.P_C_a(idx_1y);
            lhs0.const = gesol.nom_bonds.P_B_a(idx_buy) - gesol.nom_bonds.P_B_a(idx_1y);
            lhs0.slope_X = gesol.nom_bonds.P_B_b(idx_buy) - gesol.nom_bonds.P_B_c(idx_1y);
            lhs0.slope_eps = gesol.nom_bonds.P_B_c(idx_buy) - gesol.nom_bonds.P_B_c(idx_1y);

            lhst.fun_Nxz = log(gesol.nom_bonds.P_A(:,:,:,idx_sell));
            lhst.slope_Y = -gesol.nom_bonds.P_C_a(idx_sell);
            lhst.const = -gesol.nom_bonds.P_C_a(idx_sell);
            lhst.slope_X = -gesol.nom_bonds.P_B_b(idx_sell);
            lhst.slope_eps = -gesol.nom_bonds.P_B_c(idx_sell);

            % rhs = tightness
            rhs0.fun_Nxz = gesol.theta;
            rhs0.const = 0;
            rhs0.slope_Y = 0;
            rhs0.slope_X = 0;
            rhs0.slope_eps = 0;
            [a,betas_tightness(n-1),R2] = generic_nominal_reg(gesol.inflation_process, lhs0, 12, lhst, rhs0, prob_NxYz, Phi_NxYz_12mo, Y_grid);

            % rhs = f
            rhs0.fun_Nxz = gesol.f;
            [a,betas_f(n-1),R2] = generic_nominal_reg(gesol.inflation_process, lhs0, 12, lhst, rhs0, prob_NxYz, Phi_NxYz_12mo, Y_grid);

        end
    
    end

    % multiply by 100 since hpxr is not in percentages
    betas_tightness = 100*betas_tightness;
    betas_f = 100*betas_f;

end

%% old save

% prob_Nxz = squeeze(sum(prob_NxYz,3));
% 
% mean_tightness = sum(prob_Nxz(:).*gesol.theta(:));
% var_tightness = sum(prob_Nxz(:).*(gesol.theta(:).^2)) - mean_tightness^2;
% mean_f = sum(prob_Nxz(:).*gesol.f(:));
% var_f = sum(prob_Nxz(:).*(gesol.f(:).^2)) - mean_f^2;
% 
% % only compute matrix powers once (this is very slow)
% Phi_NxYz_12mo = mpower2(Phi_NxYz,12);
% 
% betas_tightness = zeros(1,4);
% betas_f = zeros(1,4);
% 
% idx_1y = find(gesol.nom_bonds.maturities==12);
% 
% [NN, Nx, NY, Nz] = size(prob_NxYz);
% temp_NxYz = zeros(NN*Nx*NY*Nz,1); % for storage
% 
% Y_NxYz = permute(repmat(Y_grid(:),[1 NN Nx Nz]),[2 3 1 4]);
% theta_NxYz = permute(repmat(gesol.theta,[1 1 1 NY]),[1 2 4 3]);
% f_NxYz = permute(repmat(gesol.f,[1 1 1 NY]),[1 2 4 3]);
% 
% lhs12 = zeros(NN*Nx*NY*Nz,1);
% lhs0 = zeros(NN*Nx*NY*Nz,1);
% 
% for n = 2:5
% 
%     idx_pnbuy = find(gesol.nom_bonds.maturities==12*n);
%     idx_pnsell = find(gesol.nom_bonds.maturities==12*(n-1));
% 
%     lhs12(:) = reshape(permute(repmat(log(gesol.nom_bonds.P_A(:,:,:,idx_pnsell)),[1 1 1 NY]),[1 2 4 3]),[],1) ...
%         - gesol.nom_bonds.P_C_a(idx_pnsell)*Y_NxYz(:);
% 
%     lhs0(:) = -reshape(permute(repmat(log(gesol.nom_bonds.P_A(:,:,:,idx_pnbuy)),[1 1 1 NY]),[1 2 4 3]),[],1) ...
%         + gesol.nom_bonds.P_C_a(idx_pnbuy)*Y_NxYz(:) ...
%         + reshape(permute(repmat(log(gesol.nom_bonds.P_A(:,:,:,idx_1y)),[1 1 1 NY]),[1 2 4 3]),[],1) ...
%         - gesol.nom_bonds.P_C_a(idx_1y)*Y_NxYz(:);
% 
%     temp_NxYz(:) = prob_NxYz(:).*lhs0(:);
%     mean_lhs0 = sum(temp_NxYz);
%     mean_lhs0_tightness = sum(theta_NxYz(:).*temp_NxYz);
%     mean_lhs0_f = sum(f_NxYz(:).*temp_NxYz);
% 
%     temp_NxYz(:) = prob_NxYz(:).*(Phi_NxYz_12mo*lhs12(:));
%     mean_lhs12 = sum(temp_NxYz);
%     mean_lhs12_tightness = sum(theta_NxYz(:).*temp_NxYz);
%     mean_lhs12_f = sum(f_NxYz(:).*temp_NxYz);
% 
%     mean_hpxr_n = mean_lhs0 + mean_lhs12;
%     mean_tightness_hpxr_n = mean_lhs0_tightness + mean_lhs12_tightness;
%     mean_f_hpxr_n = mean_lhs0_f + mean_lhs12_f;
% 
%     betas_tightness(n-1) = (mean_tightness_hpxr_n - mean_tightness*mean_hpxr_n)/var_tightness;
%     betas_f(n-1) = (mean_f_hpxr_n - mean_f*mean_hpxr_n)/var_f;
% 
% end
% 
% % multiply by 100 since hpxr is not in percentages
% betas_tightness = 100*betas_tightness;
% betas_f = 100*betas_f;